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Abstract 

The detection of primordial gravitational waves is one of the biggest challenges of the present time. The existing (Wilkinson 
Microwave Anisotropy Probe) observations are helpful on the road to this goal, and the forthcoming experiments (Planck) are 
likely to complete this mission. We show that the 5-year Wilkinson Microwave Anisotropy Probe TE data contains a hint 
of the presence of gravitational wave contribution. In terms of the parameter R, which gives the ratio of contributions from 
gravitational waves and density perturbations to the temperature quadrupole, the best-fit model produced R = 0.24. Because 
of large residual noises, the uncertainty of this determination is still large, and it easily includes the i? = hypothesis. However, 
the uncertainty will be strongly reduced in the forthcoming observations which are more sensitive. We numerically simulated 
the Planck data and concluded that the relic gravitational waves with R = 0.24 will be present at a better than 3cr level in the 
TE observational channel, and at a better than 2a level in the 'realistic' BB channel. The balloon-borne and ground-based 
observations may provide a healthy competition to Planck in some parts of the lower-i" spectrum. 
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I. INTRODUCTION 



The relic gravitational waves must have been generated by strong variable gravitational field of the very early Uni- 
verse [ij . We got used to the notion that the generation of gravitational waves requires some sort of a time-dependent 
quadrupole moment in the energy distribution of matter and fields. But, by definition, there is no such asymmetries 
in the distribution of matter and fields comprising a homogeneous isotropic universe. So, how the relic gravitational 
waves could come into existence? In the context of classical gravitational waves, the word 'amplification' is probably 
a more accurate explanation than the word 'generation'. Indeed, the pre-existing classical gravitational waves could 
be amplified in the course of time, because the nonlinear character of gravitation makes them parametrically coupled 
to the variable gravitational field of the isotropic homogeneous universe (gravitational pump field). The coupling is 
strong and effective if and when the variations of the pump field are so fast that the inverse time-scale of variations 
becomes comparable with the frequency of the wave. This is when the amplitude of the wave increases over and above 
the limits of the adiabatic law. 

Moreover, the phenomenon of relic gravitational waves is more general than this classical picture; it does not 
require the previously generated waves and does not rely on their presence. The reason is that in the quantum world 
there always exist, even in the absence of classical waves, the inevitable ground (vacuum) state fiuctuations of the 
corresponding field. One can think of these zero-point quantum oscillations as those that are being amplified, as 
a result of their interaction with the gravitational pump field. Speaking more accurately, the quantum-mechanical 
Schrodinger evolution transforms the initial no-particle (vacuum) state of gravitational waves into a multi-particle 
(strongly squeezed vacuum) state. In this sense, relic gravitational waves have been generated from 'nothing' by the 
external purnp field. This process is called the superadiabatic (parametric) amplification. For a recent review of the 
subject, see 

Under certain extra conditions, the cosmological density perturbations (dp), which combine specific perturbations of 
the gravitational field and matter fields, can also be generated by the same mechanism of superadiabatic (parametric) 
amplification. This takes place in addition to the generation of gravitational waves (gw). The theoretical ingredients 
describing the independently quantized 'tensor' {t, gravitational waves) and 'scalar' (s, density perturbations) degrees 
of freedom of the perturbed early Universe are almost identical. This results in approximately equal amplitudes of 
the long-wavelength metric perturbations representing gravitational waves on one side and density perturbations on 
the other side. Therefore, assuming that the observed anisotropics in cosmic microwave background radiation (CMB) 
are indeed caused by cosmological perturbations of quantum-mechanical origin, one expects that the contributions of 
density perturbations and gravitational waves to the lower-order CMB multipoles are of the same order of magnitude 
(see [^l and references there). 

The content of the paper is as follows. In Sec. |TT]we explain the properties of the gravitational field perturbations 
/iy (7/,x), and summarize the definitions and notations that will be used in the paper. We shall be working with the 

expansion of the field hij over spatial Fourier harmonics e**"", and with the associated Fourier coefficients Cn and 

Cn . In the rigorous quantum-mechanical version of the theory, c^„ and c^„ are the annihilation and creation operators 
acting on the quantum states. But, to simplify calculations, we are using a 'classical' version of the theory, whereby 

Cn and Cn are treated as classical random complex numbers Cn and Cn . 

In Sec.|IlT]we explain in detail how the statistical properties of the underlying cosmological perturbations, encoded 
in the probability density functions (pdf 's) of complex random variables Cn, Cn , translate into the statistical properties 
of the CMB multipole coefficients {X — T, E, B). The statistics of the CMB anisotropics, often referred to in the 
form of a loose and poetic concept of 'cosmic variance', is an important ingredient of the whole problem of deriving 
cosmological conclusions from the CMB data. We formulate the pdf 's for individual quantities af^ and for the full 
set {a^} of these quantities. The main emphasis is on the correlation coefficient pg, in the joint pdf for T and E 
components of the CMB. 

In Sec. IIVI we formulate the pdf 's for the best unbiased estimators of the auto-correlation and cross-correlation 
power spectra. Again, the main emphasis is on the estimator Dj^ of the TE power spectrum. This is because the 
gravitational waves have a distinctive signature in the region of lower ^'s: the mean value of Dj^ must be positive 
for density perturbations and negative for gravitational waves Q . The individual outcomes may have signs opposite 
to the sign of the mean value. This is why we evaluate the probabilities of finding positive (negative) outcomes for 
DJ^ in the situation when the expected value of this variable is negative (positive). 

The actual numerical amount of the TE cross-correlation at every (. is governed by the numerical value of the 
correlation coefficient pi. In Sec. |V]we derive pi for specific cosmological models with various amounts of density 
perturbations and gravitational waves. The parameters of the background cosmological model are always taken from 
the 5-year Wilkinson Microwave Anisotropy Probe (WMAP) best-fit ACDM cosmology Q, but the perturbation 
parameters are allowed to vary. The derived coefficients make it possible to build the ^-dependent confidence 
intervals surrounding the mean value of the estimator DJ^ . 
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The inherent uncertainty of the expected CMB signal, ultimately explicable in terms of the quantum-mechanical 
origin of cosmological perturbations, is exacerbated by the 'real life' effects, such as instrumental noises, foreground 
emissions, incomplete sky coverage, etc. In Sec. lVII we use the parameters of the WMAP and Planck missions to show 
how these effects broaden the size of the confidence intervals. This allows us to put to the test the existing WMAP 
TE data and make predictions for the Planck mission. 

The WMAP5 best-fit analysis suggests that the observed CMB anisotropics have no contribution from gravitational 
waves [3]. We call this conclusion the null hypothesis Hq. In Sec. IVIII we test the null hypothesis using the WMAP5 
TE data. Although the evidence for rejection of Hq is weak, within la, the data points demonstrate the tendency to 
lie below the Hq expectation curve, which is actually required by the presence of relic gravitational waves. 

Since the existence of relic gravitational waves is a necessity dictated by general relativity and quantum mechanics, 
we include gravitational waves in the TE data analysis as a compulsory ingredient (Sec. IVIIip . We consider a class 
of models which contain gravitational waves and are consistent with the available TT,EE and BB data. Although 
our physical picture is guided by gravitational waves of quantum- mechanical origin [2], the analysis is general and 
applies to gravitational waves of any origin. One has to note, though, that classical gravitational waves generated at 
later stages of cosmological evolution have wavelengths much shorter than the Hubble radius and therefore they do 
not affect the lower-order CMB multipoles. We build and study the likelihood function for R, where R is the ratio of 
gw and dp contributions to the temperature quadrupole C^^. The maximum of the likelihood function turns out to 
be at i? = 0.24. This number also determines other perturbation parameters in our class of models. Because of the 
large WMAP noises, the confidence interval for R is quite broad, so that the WMAP team's suggestion i? = is only 
IfT away from our maximum likelihood value R — 0.24. Nevertheless, we consider our best-fit model with R — 0.24 
as a benchmark model, which we believe is to be confirmed by the more sensitive Planck mission. We numerically 
simulate the future Planck data and show that the relic gravitational waves are to be detected by the TE method at 
a better than 3a level. 

The prospects of the Planck mission are discussed at some length in the last Sec. IIXI We compare the TE method 
with the more often mentioned BB method. The B-mode of polarization does not contain the contribution from 
density perturbations but is inherently weak and prone to various systematic effects. We distinguish the 'optimistic' 
and 'realistic' cases in detection of gravitational waves through the BB correlation function. It appears that the 
'realistic' BB case offers a little better than 2a detection of relic gravitational waves in the R ~ 0.24 model. We 
also show that the BB method mostly relies on the reionization era, whereas the TE correlation function is mostly 
sensitive to the era of recombination. 

The Conclusions summarize our findings and emphasize the necessity of concentrated efforts in the race for discovery 
of relic gravitational waves. 



II. GRAVITATIONAL FIELD PERTURBATIONS 



Let us recall that the gravitational field of a slightly perturbed flat Friedmann-Lemaitre- Robertson- Walker (FLRW) 
universe is described by 

ds^ = -c^dt^ + a^{t){Sij + h,j)dx'dx^ = 0^(77) [-^7^2 (^g^^ ^ h,j)dx'dx^], 

where t-time and 7y-time are related by cdt = adrj. The metric perturbation field hij{ri,x) can be expanded over 
spatial Fourier harmonics e**"": 



P,j (n) tin (?7)e'"-^ c„ + p* (n) K {v)e'"^"" cl 



(1) 



where n is the dimensionless time-independent wavevector. The dimensionless wavenumber n is n = {Sijn^n^y^^ . 

The dimensionless n is related with the dimensionful k hy k = n/2///, where Ih = c/ Hq is today's Hubble radius. 
The value nn = 47r labels the wave whose today's wavelength is equal to the 'half of the size of the Universe', that 
is, equal to today's Hubble radius. The use of n, rather than fc, is very convenient as it guides the evaluation of the 
CMB multipoles: cosmological perturbations with wavenumber n are mostly responsible for the CMB temperature 
anisotropics at the multipole £ ^ n. 

The general form of Eq. ([1]) is common for metric perturbations hij(ri,x) characterizing gravitational waves and 

s 

density perturbations, but gw and dp differ in the explicit content of the two polarization tensors Pij (n) (s = 1,2). 
In the gw case the field hij consists of two 'transverse-traceless' components, whereas in the dp case it consists of the 

'scalar' and 'longitudinal- longitudinal' components (see, for example, Q). The gw mode functions hn (v) satisfy the 
perturbed Einstein equations linearized around an FLRW background. In the dp case the gravitational mode functions 
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hn (^) are necessarily accompanied by matter mode functions. For many models of matter, the full set of dynamical 
equations for density perturbations leaves independent only one out of two polarization components s = 1, 2. 

In the rigorous quantum-mechanical version of the theory, the quantities hij are quantum-mechanical operators 

acting on quantum states. The and are the annihilation and creation operators, satisfying the commutation 

relations [c„, cji] — (5s's(5('^)(n— n'). In the case of density perturbations, the same operators and are inherited, 
through the linearized Einstein equations, by the matter fields. According to the physical formulation of the problem, 
the initial quantum state of the quantized perturbations is the ground state |0) of the corresponding Hamiltonian. The 

ground state satisfies the condition |0) = 0. For gravitational waves, the normalization constant C is C = VTGttZpi, 

whereas for density perturbations C = \/2AttIpi, where Ipi = (^Gh/c^) is the Planck length. Obviously, if the Planck 
constant h is artificially sent to zero, the initial zero-point quantum fluctuations, as well as the field hij itself, vanish. 
To simplify calculations, we will be using a 'classical' version of the theory, whereby the quantum-mechanical 

operators c„ and cj are treated as classical random complex numbers c„ and Cn. We are loosing some delicate 



s * 



quantum aspects of the rigorous theory but they are not our main concern. The statistical properties of €„, c 
define the statistical properties of the classical field hij ^ and, through the linear radiative transfer equations, they 
determine the statistical properties of the CMB anisotropics. 

Since the Schrodinger evolution of the initial vacuum state into a squeezed vacuum state retains the Gaussianity of 
the state (while developing strongly unequal variances in the amplitude and the phase), we are justified in assuming 
that each individual complex coefficient is described by a complex pdf with a zero mean and unit variance. 
If necessary, one can imagine that the space of wavevectors n is discretised into a fine 3-lattice. Moreover, since the 
different modes n and polarization states s are not supposed to 'know about each other', we are justified in assuming 
that the joint distribution for all coefficients {c^} (that is, for the full set including all possible values of n and s) is 
a complex multivariate Gaussian pdf comprising a product of individual Gaussian pdf's [5, 6]. 

The first and the second moments calculated with the joint pdf amount to 

(c„) = (c„*) = 0, {CncJ) = 6,s'S^^\n - n'), (c„c„,) = {KcJ) = 0. (2) 

Since the coefficients c„ obey a joint multivariate Gaussian (normal) distribution, the higher order averages are 
expressible in terms of the second moments. 

It is clear from Eq. ^ that the mean value of the field hij is zero, but the variance is non-zero and is calculated to 

be 

OO 

^{h,,{r^,^)h'^rj.,^)) ^ I ^h'{n,r,). 







The function 



h'{n,v)^f^n^ E (3) 



s=1.2 



gives the mean-square value of the gravitational field perturbations in a logarithmic interval of the wavenumbers 
n. It is called the metric power spectrum. The spectrum ([3]) refers to gravitational waves or density perturbations 
depending on whether this quantity is calculated from gw or dp mode functions. 

The CMB calculations are usually being done numerically, whereby the initial conditions for differential equations 
are imposed at times deep in the radiation-dominated era 0, Q. The wavelengths of relevant modes n are much 
longer than the Hubble radius at that time, and the primordial power spectrum /i^(n) is a smooth function of n. The 
early pump fields (scale factors) a(ry) which are power-law functions of rj generate primordial spectra h'^{n) which are 
power-law functions of n : 

h^n) (gw) = Bln^\ h^{n) {dp) = B^n"""!. (4) 

The quantum-mechanical evolution of the initial vacuum states of gw and dp fields results 2] in approximately equal 
primordial amplitudes B^, B^, and also in the relationship rit ~ Ug — 1. Although, in general, the four numbers 
Bf,B^,nt, Us specified at the radiation-dominated stage do not contain enough information to impose all the starting 
conditions for numerical calculation of the ensuing GMB anisotropics, some additional simplifying assumptions 
[1| allow one to do that. 
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At the end of this section it is necessary to warn the reader that according to the claims of inflationary theory 
the amplitudes of density perturbations should be arbitrarily large in the limit of models with small values of the 
early Universe parameter e, where e = —H/H^. Specifically, the incorrect (inflationary) treatment of quantized scalar 
perturbations has led to the proposition that the power spectrum of scalar metric perturbations (often denoted by ^ 
or TZ) should be divergent as 1/e, right from the very early time when initial conditions in the high-frequency regime 
of the perturbations were imposed. Inflationary theory starts its line of argument with vacuum fluctuations of the 
scalar field (inflaton) in de Sitter space time and finishes with infinitely large scalar metric perturbations in the same 
de Sitter spacetime. This predicted divergency of density perturbations in the limit of e = {rig = 1, nt = 0) is 
customarily being substituted by the claim that it is the amount of primordial gravitational waves that is expected 
to be very small. This is expressed in the form of small values of the 'tensor-to-scalar' ratio r: 

r = 16e = ~8nt. (5) 

In particular, according to the inflationary Eq. ([5]), the interval of a 'standard' (de Sitter) inflation, i.e. e = 0, 
generates a zero amount of primordial gravitational waves, i.e. r = 0. This is in full contradiction with the mechanism 
of superadiabatic (parametric) amphfication 0. Certainly, we are not using Eq. ®, neither as a valid theoretical 
result nor as an element of CMB data analysis. For a detailed critical discussion of inflationary predictions and the 
data analysis based on inflationary theory, see [3] . 

III. ANISOTROPIES IN THE COSMIC MICROWAVE BACKGROUND RADIATION 
A. Characterization of the radiation field 

The radiation field is usually characterized by four Stokes parameters (/, Q, U, V), where / is the total intensity of 
radiation, Q and U describe the magnitude and direction of linear polarization, and V is the circular polarization. 
The Stokes parameters can be viewed as functions of (t, x*, j/, e*), where v is the photon's frequency, and e* is a unit 
vector in the direction of observation (opposite to the photon's propagation). In a given space-time point (i, a;*) and 
for a fixed v, the Stokes parameters are functions of 6, cf), where 0, (f> are coordinates on a unit sphere indicating the 
direction of observation: da^ = gabdx°'dx^ = dO^ + s\t? 9d4>^ . 

The Stokes parameters of the radiation are components of the polarization tensor Pab 9] which can be written as 

I + Q -{U-iV)iiine \ 
-{U + iV)iiine (/ - Q) sin^ 6* j ■ 

Under arbitrary transformation of coordinates (0, 0), the components of Pab{0, (f) transform as components of a tensor, 
but some combinations oiPab{0,4>) and its derivatives remain invariant. 
Two algebraic invariants are given by 

/(0, 0) = g''\e, ^)Pab{e, cl>), V{9, 0) = ze'^'iO, ^)Pab{e, cl>), 

where e°'^{6^ (f) is a completely antisymmetric pseudotensor 

ab^f -(sin0)-i\ 
^ (sin^)"' ) ■ 

Now one can single out the symmetric trace-free (STF) part P^b^ of Pab'- Pab{0, 4>) — \lgab — ^Vcab + Pab'^ ^ 

pSTF^lf Q -Usin0 \ 
2 V -C^sin6' -Qsin^e J ' 

Two other invariants require the use of second covariant derivatives of P^b^^ ■ It can be shown that there exist only 
two linearly independent invariants that can be built from second derivatives 0: 

E{0, Ct>) = -2(pSTF);a;6^ ^(^^ ^) ^ ^^(^pSTFyb;d^a^ 

In the above expression ";" denotes covariant differentiation on the sphere. The quantities / and E are scalars, while 
V and B are pseudoscalars. 
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In CMB applications one usually ignores V, as well as the unperturbed part of /, and regards the remaining 
quantities as measured in the temperature units, specifically in micro-Kelvin (/iK). These quantities (temperature 
and polarization anisotropies) can be expanded over ordinary spherical harmonics YimiS, 0): 



no, 



(6a) 



=0 m=-i 



EE 

e=2 m=-l 



=2 m=-e 



(£ + 2) 



(6b) 



The ^-dependent factors in front of af and in (I6bp have been introduced in order to work with definitions 
consistent with previous literatures [1], p^ . 

Since /, E, B are real functions on the sphere, the complex multipole coefficients in ((6a|l and (|6bp satisfy the 
following conditions 



{X^T,E,B, ~£<m<£). 



(7) 



The set of multipole coefficients {ojjyi, af^^, af^^) completely characterizes the CMB field that we are studying. 



B. Statistical properties of the CMB field 

The radiative transfer equations relate the set of random coefficients with the multipole coefficients 



(« 



T 



). As a consequence, the multipole coefficients {X 



T, E, B) inherit the randomness of €„. 



The calculations usually begin from a single Fourier mode of perturbations considered in a special reference frame 
whose z-axis is parallel to the wavevector n. This is how one arrives at the multipole coefficients denoted a^(n, s). The 
quantities afj^{n, s) depend in a complicated but calculable manner on the metric (plus matter, in the case of de nsity 
erturbations) mode functions and the background cosmological model [H [H, [H, H, [S [H, E E S [U, [11, 
!23j . The next step is the generalization of this result to the set of all Fourier components with arbitrarily directed 
n's. This requires the use of the Wigner rotation coefficients D'^^,{h), where n = n/n. The Wigner coefficients 
describe the rotations of individual special frames associated with n's to one unique frame of the observer. The 
Wigner coefficients satisfy the orthogonality relation [23| 



i2£+l) / Dl^, {n)D^M,{n) d ilf, ^ AttSilS. 



JmMOm'M' 



(8) 



A detailed analysis shows fs^ that the resulting multipole coefficients af^ are related to the random coefficients €„, 
both for gw and dp, as follows 



C 



(27r)3/2 



d^ 



n 



E 



J2 [alAn,s)Di^,{h) c„ +(-l)™af„:,(n, s)i? 



,(n) 



(9) 



The linear character of this relation implies that the probability density functions of individual multipole coefficients 
af^ are zero-mean Gaussian pdf 's, and the joint distribution of all multipole coefficients is described by a multivariate 
Gaussian pdf. The effects of nonlinearities and quantum loops of the underlying metric field are very small, and we 
ignore them. The nonlinear effects in the radiative transfer equations are possible, but apparently they are also small. 
The individual distributions are marginal distributions derivable from the joint pdf. Before formulating these pdf's 
explicitely, we shall use Eq. ^ and Eq. ^ in order to find second moments of the distributions. 



1. Auto-correlation and cross-correlation functions 

It follows from Eq. ^ and Eq. ^ that all linear averages of vanish. 
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The quadratic averages {afjy^a-^^,} do not vanish. They are called the correlation functions. One gets an auto- 
correlation function ii X = X', and a cross-correlation function ii X X' . 

Let us start from the auto-correlation functions [H, [13, El, H^l , 3, 0- Using Eqs. ([21 El El) one can show that 

{af^af,*^,) = Smm'Su'Cf^ , (af^a^^,) = {af^af,^,) =0, {X ^ T,E,B), (10) 

where the XX power spectrum Cf^ is given by 0] 

^^"^ = 2.^(2£+l) /"^" E E \'^Lin.s)\\ (11) 



= 1,2, 



Both, density perturbations and gravitational waves, produce distinctive CT^ , power spectra, but the Cf^ is 

zero for density perturbations and Cf^ ^ for gravitational waves [1, [ToT [25l . [26| . Specifically, all the coefficients 
vanish in the dp case. In other words, density perturbations do not produce the _B-mode of polarization 



ialli, 3- 

Let us now turn to the cross-correlation functions. In general, there exists three cross-correlation functions, namely 
TE, TB and EE. However, the TB and EE functions vanish identically in the dp case, as density perturbations 
do not produce the i3-mode. In the gw case, the left and right circular polarizations of gravitational waves give 
contributions of opposite signs into the TB and EB functions. It is very natural to expect that the relic gravitational 
waves with both polarizations had been generated in equal amounts. Then, one concludes that in the gw case the TB 
and EB functions do also vanish 3- (However, there exist parity-violating processes in the later Universe, see 
for example [11] ■) Therefore, we shall only deal with the TE cross-correlation function. 

Using Eqs. (0 11|51), we arrive at 



(12) 



where the Cf^ is given by 

'^^^ " 47^2(2^+1) / E E [«L(":S)afr^('^>s) +«rm('i>s)aL(":s)] . (13) 

s— 1,2 m——i 

Note that one can derive Eqs. (fTOl \12^ as soon as the averages ([2]) are given, independently on whether the actual 
distributions are Gaussian or not. Also, it is seen from Eqs. (|10[ [T^ that only the products of coefhcients with £ = i' 
and m = m' are important. 

Unlike the auto-correlation power spectra Cf^ which are strictly positive functions, the values of the TE 'power 
spectrum' Cf^ can be positive or negative. Most importantly for our further discussion, in the region of lower ^'s 
the function Cf^ is positive for density perturbations and negative for gravitational waves Q. It can also be shown, 
using Eqs. (HHni), that 

{crr < crcr, (m) 

which is the analog of the Cauchy-Bunyakovsky-Schwarz inequality. 

Since density perturbations and gravitational waves are independent fields, they contribute to the total CMB power 
spectra additively: 

^' = ^' (dp) + ^' {gw}. (15) 
2. Probability density functions for individual multipole coefficients 

The Gaussian nature of Cn, together with Eq. ([51), translate into Gaussian pdf's for individual multipole coefficients 

'^em- 
It is convenient to start from distributions for real and imaginary parts of a^, 

X _ X{r) . X{i) 

It follows from Eq. ^ that 

„-fW _ n X(r) _ / -.Nm X(r) , _ / -|^m+1 
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Thus, the set of multipole coefRcients af^, where i > m > —£, is fully equivalent to the set of real quantities 

afj'''^ , afj[^ , afj^^^ , where £ > m > 1. We will be using index c to denote all quantities as afjf'' where c = r or i if 
m > 1, and c — r if m = 0. 

Each individual coefficient afj^'^^ satisfies a zero-mean normal distribution 

J y^ern ) - r— X(c) ^ 

The variances, i.e. the squares of standard deviations crfj^'^\ are the second moments already analyzed in the previous 
subsection. Comparing with Eq. (|10p one finds 

(rr^(^)\2 _ , X(r) X(r\ _ ^XX 
X(r) X(i 



X^-fm "fori I — 



Denoting af = ^Gf-^, one can write: cr|^ = af and u^}^^ = crf^^^ — uf In other words, the standard 

deviations of the (21 + 1) random variables , af^j^'\ '^fm \ where £ > m > 1, are all equal and amount to one and 
the same number /V2. 

The equivalent statements follow from the pdf for complex multipole coefficients af^ (m > 1). Each individual 
obeys a complex zero- mean normal distribution [5, 6] with the standard deviation erf: 



1 



fiafj = — ^e-('^-'^-)/(-^'r. (16) 

7r{(Tf )2 



3. Joint probability density function for multipole coefficients 

Although the joint pdf for the set of all coefficients Cn is a factorizable Gaussian distribution, the radiative transfer 
equations make similar factorization impossible for some subsets of all multipole coefficients af^. Specifically, aj^ 
and afi^' 'do know about each other', when £ — £' , m — m' . Their joint pdf is not factorizable and their correlations 
do not vanish, as was made clear in Eq. (fT2ll . 

Let us denote sets of coefficients by curly brackets. The full set of all multipole coefficients {a^} is 

{af„,} = {afj X^T,E,B; £^2,3,---, -£,-■■,£). 

In this set we did not include the potentially nonzero coefficients o^q, aj^^, and hence, from now on, we restrict our 
analysis to multipoles £ > 2. 

Each member of the set {a^} is a linear combination ^ of complex coefficients whose joint distribution is 
a zero- mean normal pdf. This guarantees that the joint distribution for the set {ct^} is a zero-mean multivariate 
normal distribution [5|. Such distributions are completely characterized by their second moments. For the physical 
reasons explained in Sec. lIIIB II we assumed that there is no TB and EE second moments, so that the B components 
'do not know' about T and E components. As for the non- vanishing second moments, they were calculated in Eq. (jlOp 
and Eq. (HH). 

Combining the above statements into a formula, we can write the joint pdf for the full set {afj^}'- 

f {{^Jrm '^frm ^fm}) — f {{'^Jrm ^fm}) f ii'^fm}) ■ 

Here, 

oo e 

/({aL})=n n 

£=2 m=-t 

where each of the functions f{af^) (for fixed £ and m) is a single-variable normal distribution ([16]) with X ^ B. As 
for the T,E part of the joint pdf, we are interested in the subset of all a.J^,af,^,, in which £ = £\m — to'. This 
part of the joint pdf is the product of the functions /(a^, o,fm)i for each fixed pair of £ and to, which we are set to 
formulate. 
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The joint distribution for complex a^, af^ {m 7^ 0) is a bivariate normal distribution which can be written as 



m ' 



exp 



T T* E E* ( T E* j_ T* E \ 



T^2 



E\1 



(17) 



In the case m = 0, the joint pdf reduces to a bivariate normal distribution for real quantities a 



1 



: exp 



1 



-.T „E-\ 



jam) , (afo) _ ^P^a^a 



(18) 



The aj and af in Eq. (fT7|l and Eq. (|18p are the standard deviations of aj^ and af^, respectively. The pe in 



Eq. (fT7|) is the correlation coefficient, defined by 

L(„T fjE* , fjT* E \ 



Pi 



r^TTr>EE 



(19) 



Like the standard deviation af, which, for a given is one and the same for all m's, the correlation coefficient pg 
is also one and the same for all m's, i > m > ~£. For a bivariate normal distribution the correlation coefficient 
obeys the condition \pi\ < 1, which essentially is Eq. p4|) . Concrete numerical values of p^'s depend on the type of 
underlying perturbations and cosmological background model. 

The complex pdf ITTl) with m ^ is equivalent to the product of joint pdf's for real and imaginary parts of 

(a^, af^), i.e. for {cij^\afj^^) and {aJ^nKaf^^) with £ > m> 1 [5i]. Introducing the label c = r or i one can write 
the common formula 



1 



Tiajafyjl-pj 



exp ■ 



1 



i^-pi) 



Tic) E{c) 



{of? 



(20) 



The quantities {aj , erf , af , pg} participating in and defining the pdf's are related with the CMB power spectra 
calculated in Sec. IIIIBI 



/^T\2 r'TT („E\2 



EE 



We ) = '-'I 



BB 



TE 



TE 



(21) 



We will use these relations in our further discussion. 



IV. ESTIMATORS OF THE CMB POWER SPECTRA AND THEIR PROBABILTY DENSITY FUNC- 
TIONS 



The quantum- mechanical origin of cosmological perturbations, which we represented by randomness of c„, does not 
allow one to predict a unique CMB map, i.e. one concrete set of multipole coefficients af^- The most that theory 
allows is to formulate probability density functions, and we have done this. The pdf's are fully characterized by the 
set of numbers {aj , af , af , p^} or, equivalently, by the power spectra (|2ip . Theory can predict these parameters, 
while observations can check whether these parameters are indeed as predicted. 

To find a parameter from observations, one normally uses an appropriate random variable (statistic) called an 
estimator. The estimator is unbiased if its mean value is equal to the parameter, and is called the best if its variance 
is as small as possible. A concrete value of the estimator calculated on the actually observed data (in our case, on 
the set of actually observed multipole coefficients a^) gives an estimate of the parameter. We will denote by symbol 



-D^^ the best unbiased estimators of the power spectra Cf^ and will consider separately the auto-correlation, 
X = X' and cross-correlation, XX' = TE, estimators. 



A. Estimators of the auto-correlation power spectra 



It is clear from Eq. (fTU]) that there exists plenty of unbiased estimators of Cf''^ . Each of 2£ -|- 1 estimators (random 



variables) (a 



with £ > m > —£ calculated on the observed multipoles provides an unbiased estimate of C 



XX 
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However, there exists only one best unbiased estimator. It is given by 

^ ' X X* 

= 1 X! {"'emO'frn) ■ 

m——i 

The expectation value of this estimator is, as required, 

(Dr) = ^ E («f™«f™) = i-fr - cr- (22) 



The standard deviation is ADf^ = J{{Df^Y) - {Df^^, amounts to [2 



ADf^ = (af)Y^. (23) 

This is the smaUest, but still a non-zero, variance among all possible quadratic estimators [30]. (Formula provides 
one of interpretations of the term 'cosmic variance'.) 

We have access to only one realization of the inherently random CMB field, that is, we can deal with only one 
observed set of multipole coefficients a^. Under certain conditions, the access to only one realization of the stochastic 
process is sufficient for extraction the true parameters of the process, with probability arbitrarily close to 1, from 
the data. When this is possible the stochastic process is called ergodic. However, there is no ergodic processes on a 
2-sphere (our sky) [3^. We will always be facing some uncertainty surrounding the parameters extracted from the 
observed CMB field. 

The variance of an estimator at each i is just a number, and is not sufficient for proper handling of the data. The full 
available information is contained in the pdf of the estimator. In order to derive the pdf of Df-^ , {XX = TT, EE, BB) , 
we will first rewrite it as 



m— 1 



It is clear that the estimator Df^ is just the sum of squares of (2£ + 1) independent normal variables. Thus, for 
every fixed £, the estimator obeys a distribution with 2£ + 1 degrees of freedom 

In what follows. Sec. IVI Bl we will be dealing with the problem of observations on a cut sky. The cut sky limitation 
effectively means that, for a given £, one has access not to all 2i+l degr ees of freedom, but only to a smaller number 
n, n = {2£ + l)/sky, where /sky is the observed fraction of the sky [3l|, [H, H^l- In preparation to this problem, we 
are writing down the pdf for Df^ with n degrees of freedom, Df-^{n), where n is not necessarily equal to 21+1. 
Denoting V = nDf-^{n)/{afY, one writes Q 

„T/n/2-l -y/2 

^ 2-/.r(„/2)(.,v). . P") 

where T is the Gamma-function. Certainly, when n = 2£ + 1, the mean value and the standard deviation of the 
variable Df^ derivable from the pdf ((24l) . coincide with the already written expressions (l22t and ((23l) . 



B. Estimators of the cross-correlation TE power spectrum 



The TE cross-correlation is in the focus of this paper, so we shall explore the TE estimator in depth, without 
sparing details. We know that the Cf^ correlation function should be positive, if caused by density perturbations, 
from i = 2 and wp to £ 53; and it should be negative, if caused b y g ravitational waves, from £ = 2 and up to 
£ fa 150 0, a further study of the TE correlation was conducted in [34| . But Cf^ is only the mean value of the 
correlations. In a particular realization of the random field, i.e. for the actually observed af^, the Dj^ outcomes 
can be negative, despite the fact that the theoretical Cj^ is positive, and vice versa. Suppose negative correlations 
are observed in a set of multipoles. Is it necessarily a signature of gravitational waves, or it may be a statistical fiuke 
of density perturbations alone? Suppose positive correlations are observed in a set of multipoles. Are they 'positive 
enough' to be explained by density perturbations alone, or they require an admixture of gravitational waves? This is 
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the kind of fundamental questions of the TE approach that we need to answer [2] , and they can only be answered on 
the basis of the probability density functions for the estimators. 

It is clear from the previous discussion and Eq. ([T2|) that the best unbiased estimator Dj^ is given by 



nTE _ 1 / T E* . T* E 



It can also be written as 



^r = ^f-l? + E(-t!+-^a)) (25) 

\ m=l / 

where 

As before, we denote quantities ([26]) by u^^, where c = r or i, m > 0. Obviously, as the coefficients a^^, are real, 

(i) 
«£0 



Each of the quantities vf^ is an unbiased (but not the best) estimator of Cj^ , as follows from the averages based 
on the joint pdf pi)) : 



\"hnl - P^^l - ^£ ' ^^£m - V + ^^^^ '^^ ^ ^ + l^f i ■ 

In order to analyze the pdf of the estimator Dj^ in it is instructive to study in the pdf's for quantities w^.^^. In 
Appendix [X] we present a detailed discussion of these pdfs and some of their relevant properties. 

Based on the above analysis and the results from Appendix |^ we are now in the position to derive the pdf for the 
best unbiased estimator Dj^ . Firstly, the mean value and the standard deviation of Dj^ are known from the joint 
TE pdf: 



1 ^ 

^ I T E* . T* E \ T E r<TE (r,r,\ 

rpYy 2^ (a^mflfm + ain^ag^) ^ pi>ag cr<> = , {11) 



2(2^+1) 



Aor ^ ^{{Drr) - {Dry = -^-"y |^ = f' \ti = ^fj- ^^'^ 

Clearly, ADj^ is smaller than Aw^^ by a factor \^2£ + 1. The minimum of ADj^ is sX pi = Q and the maximum 
is at Ip^I = 1. The maximum is only V2 times greater than the minimum. For growing i and fixed crj, af ^ pi, the 
standard deviation is decreasing as l/^/2£ + 1. 

To derive the pdf f{Dj^) one notes (see Eq. (|25|) that the variable Dj^ is essentially the sum of 2^ + 1 products 
of normally distributed variables from two classes, T and E, with standard deviations aj and af, respectively. Since 
we will be interested in the situation where the number n of accessible degrees of freedom is smaller than 2£ + 1, we 
have to make the derivation more general. It is instructive to start from a somewhat abstract problem. 

Let the three variables x,y, z be defined by 

x = Xf + Xi + --- + Xl y = Y,^ + Yi + --- + Yl z = X^Y^ + X^Y^ + • • • + X„y„, (29) 

where Xi and Yi [\ < i < n) are real variables. We assume that within each class, X and F, the variables are 
statistically independent and each Xi and Yi obeys a zero-mean normal distribution with standard deviation a-^ 
and CT^, respectively. We also assume that each pair Xi^Yi obeys a zero-mean bivariate normal distribution with 
correlation coefficient p. Then, the joint pdf f{x,y, z) is a Wishart distribution with n degrees of freedom Q: 

ff \ I 1 {xy-z^y 

f(x,y,z) - 



4(i-p2)(^x^y)2y ^i/2r(|)r(ii^) 

/ I / X y 2pz 
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From this distribution, by integrating over variables x and y, one derives the pdf f{z): 



/ 1 \ p z \ ( \z\ \ 

^^'^ " Y 2«-i7r(a^a^)"+i(l - P^) r(f ) \ 1 - a^a^ ] [{1 - p^)cjXa'^ ) ' 

where Kr^-i_ is the (^i^)-order modified Bessel function. In this derivation we have used (2.382.2) and (3.471.9) from 
Rcf. [35]. In the special case n = 1, and denoting the pair Xi,Yi by ajQ^af^ or \/2aJ^\ \^af^\ (for fixed m > 1, 

c = r or i), one finds that the pdf (|30|) returns to the pdf (jAl[) for each w^^. 

The pdf for the estimator Dj^ (n) with n degrees of freedom, where n is not necessarily equal to 2£+l, follows now 
from Eq. (|5D|) . Comparing Eqs. ([^ we write Dj^{n) = z/n and, then, the pdf for Dj^{n) (omitting argument 
n): 



DT^\\~ f n \— I 1 1 



TE\ _ 



2 ; w^fj \<^-pi)m 



Certainly, Eqs. ([271128]) are consequences of this pdf, when n — 2i+l. In Appendix [Blwe discuss some useful properties 
of this pdf. These properties are used in our consequent analysis, in particular for determining the confidence intervals 
of DJ^ for specific cosmological models. 



V. CORRELATION COEFFICIENT pe AND ESTIMATOR Dj IN SPECIFIC MODELS 

The sign and value of the TE correlations, as well as their statistical distribution, crucially depend on the sign 
and value of the correlation coefficient pf. As we know, at lower £'s, pi should be negative for gravitational waves 
and positive for density perturbations. However, specific numerical values of pi are different in different models. It is 
difficult to find pi analytically, so we are using numerical calculations to find it numerically. 



A. Numerical evaluation of the correlation coefficient pi 

Assuming that the correlation functions Cj^ , Cf^, Cj^ are known from numerical modelling, the pi can be found 
from the defining expressions (fT9)) . (fT5)) : 



^crcf^ ^/{cridp) + crigw)) {cf^idp) + cf^^gw)) ■ 

In our numerical calculations we work with a single cosmological background model, which is the 'best-fit' ACDM 
cosmology based on the 5-year WMAP observations and other data sets . The parameters of the background model 
are 

r^fo/i^ = 0.02265 ±0.00059, f^c/i^ 0.1143 ± 0.0034, , . 

r^A = 0.721 ±0.015, T^ejon = 0.084 ± 0.016. ^ ' 

The associated value of the Hubble parameter is /i = 0.701 ± 0.013. In calculations, we are using the central values of 
these parameters. 

As for the perturbations, the WMAP5 'best-fit' analysis suggests no contribution from gravitational waves, whereas 
density perturbations are characterized by the spectral index Us and the amplitude [J, [sl] : 

Us = 0.960l:g:gi3, ^ (2-457tg:E]|^3) X 10"^ (at k = 0.002Mpc~^) . (34) 

The WMAP5 'best-fit' parameters ([33]), ([53]) imply the temperature quadrupole 



27r 



1160 {pK^). (35) 
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FIG. 1: The correlation coefRcient pi as a function of I in models with different amounts of gravitational waves. The solid line 
denotes the best-fit ACDM model with density perturbation alone, the dashed line is for the R = 1 mixture of dp and gw, and 
the dotted line is for gravitational waves alone. 



Being guided by the quantum-mechanical theory of cosmological perturbations (see Sec.|T|, we include gravitational 
waves in our analysis from the very beginning, as a necessary ingredient of CMB anisotropics at lower £'s. We allow 
the dp parameters ([M]) to vary and strive to achieve the total theoretical TT, EE, TE spectra as close as possible 
to the observed data. Gravitational waves are parameterized by the primordial spectral index rit, where nt = — 1 
(unless some other rit is chosen for purposes of illustration), and the amplitude parameter i?, where R is the ratio of 
the gw and dp contributions to the temperature quadrupole: 



To avoid any ambiguities in definitions and derived upper limits, we do irot use the 'tensor-to-scalar ratio' parameter 
r (see |2|). For very rough comparison of results, one can keep in mind the approximate relationship r ~ 2R, which 
follows from the study |^]. In power spectra calculations we rely on the CMBFast code j38j (Alternatively one can 
use the CAMB which gives consistent results). The correlation coefficient pg is then found from Eq. p2|) . 

In the way of illustration, we plot in Fig. [T] the pi for three representative sets of perturbations on the background of 
ACDM cosmology (|33p . The solid curve shows pg for density perturbations alone with best-fit parameters ([M)) . The 
dotted curve shows pi when gravitational waves alone are present with nt = and the amplitude normalization such 
that the temperature quadrupole is equal to the best-fit WMAP5 prediction ([55]) . The dashed curve shows pg arising 
in the R — 1 mixture of density perturbations (rt^ = 0.96) and gravitational waves {nt = 0) with total quadrupole 
equal to the WMAP5 best-fit value ([55)1 . 

As expected Q and seen in Fig. [U pg is positive for dp and negative for gw at lower multipoles, £ < 50. In the 
intermediate case R = 1, pg is small \pe\ < 0.2, as the dp and gw contributions largely cancel each other. At higher 
multipoles, £ > 100, the R = I curve approaches the dp curve, as the role of gravitational waves diminishes in 
comparison with density perturbations [3, 40]. The oscillatory features at very low multipoles £ ~ 10 — 20 (clearly 
seen in the solid curve) are real. They are a reflection of effects of Q\ in combination with the reionization era 
characterized by the optical depth Treion- 

Knowing pi as a function of £ one can build and analyse pdf 's for the Dj^ estimator in specific models. As we 
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Multipole: I 



FIG. 2: The {£ + l)Cf^ /2it power spectrum in units of (/iK^) for the best-fit WMAP5 model with no gravitational waves. 



show below, the sign difference in pis for density perturbations and gravitational waves is crucial for telling one sort 
of perturbations from another. 



B. Probability of negative values of the estimator Dj 



In the center of our attention is the WMAP5 best-fit model ([55]) . ([M)) with no gravitational waves, R = 0. The 
mean value of the estimator Dj^, i.e. the TE power spectrum, is shown in Fig. O To assess the consistency of this 
graph with the actually observed data points 0, |4l| (see also [i^), many of which are negative at £ < 50, we start 
from the pdf for this model. We first ignore noises and systematic effects, but will take them into account in Sec. IVII 

Our derivation is based on Eq. (pij) . with n — 2£ + 1, and pi, as plotted by solid line in Fig. [TJ In Fig. Owe show 
f{Dj^) for a series of multipoles, ( — 2, 10, 20, 30, 50, 100. When £ increases, maxima of individual pdf 's move from 
positive to negative values, crossing zero at £ « 53, in agreement with Fig. [2] The spread of individual pdf's becomes 
narrower with increasing £ (note increasing horizontal scale on different panels), in agreement with the decreasing 
standard deviation ADf^ cx afaf /y/{2£+ 1), as reflected in (p8| 
approximation to the exact pdf's, 



We also plot by dashed curves the Gaussian 



TE\ 



1 



/2^ADT^ 



exp 



TE 



2{ADJ^Y 



(37) 



As could be anticipated, the Gaussian approximation becomes progressively better with the increasing £. 

We shall now draw various confidence intervals surrounding the mean value at each £. The confidence regions are 
determined by the shortest interval between the upper {Dj^)u and the lower {DJ^)l boundaries enclosing a given 
surface area under the pdf, as discussed in Appendix VK\ 

In Fig. m we draw confidence intervals for all three models considered in Fig. [TJ Namely, the best-fit model with 
density perturbations alone (left panel), the R = 1 mixture of density perturbations and gravitational waves (middle 
panel), and a model with gravitational waves alone (right panel). It is seen on the left panel that the 99.7% confidence 
region is above the zero line for 15 < ^ < 29. The 95.4% confidence region is above the zero line for £ < 37, while 
the 68.3% confidence region lies in the positive territory for all £ < 45. In contrast, the right panel for gravitational 
waves alone shows that all plotted confidence regions are well below the zero line for ^ > 25. 

Finally, in Fig. O we show the probability P{Dj^ < 0) of finding negative data points in models with varying 
amounts of density perturbations and gravitational waves. In all combinations of dp and gw^ the spectral indices are 
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FIG. 3: Probability density functions f{Dj'^) at £ = 2, 10, 20, 30, 50, 100 for the WMAP5 best-fit model JSl)- The dashed 
curves show the Gaussian approximation to the exact pdf's. 




FIG. 4: The 68.3%, 95.4%, 99.7% confidence intervals of the variable DJ^ for three representative models shown in Fig. [T] In 
all three panels, the solid black curves are the mean values of the statistic {£+l)Df^ /2tv, i.e. the power spectra {£+l)Cf^ /27v. 
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FIG. 5: The probability of finding negative TE data points in different models. The curves are ordered from the bottom to 
the top by increasing R, R = 0.0, 0.05, 0.1, 0.3, 0.5, 1.0 with gw alone at the very top. 

Us = 0.96, Tit = and the total quadrupole is equal to the best-fit value ((35| . The black solid curve is for the WMAP5 
best-fit model vifith no gravitational waves, R = 0. It is seen from the graph that in the case of density perturbations 
alone, the probability P{Dj^ < 0) is pretty small, especially in the range of multipoles 10 < ^ < 30. The lowest 
point is at = 20, with P{DJ3^q < 0) = 9.33 x 10^^. Obviously, when the amount of gravitational waves increases, 
the P{Df^ < 0) approaches 1. 

VI. NOISES AND OTHER SOURCES OF INCREASED UNCERTAINTY 

The spread of the probability density functions discussed so far is the unavoidable consequence of the inherent 
randomness of the sought after signal a^. Additional randomness, and additional uncertainty in the parameter 
estimation, is brought in by the process of measurement. The output (observed) readings are contaminated by 
instrumental noises and various foreground radiations, such as the synchrotron radiation, thermal emission from dust, 
unresolved extragalactic sources, etc. [1^ IH, IH, |4^, IH, 113, 13 . It is common to treat the residual systematic causes 
as 'effective noises' alongside with the 'genuine' instrumental noises. In addition, the incomplete sky coverage and the 
sky map cuts [sil . [3^ . [ssl |49 !| introduce extra loss of information in comparison to what is available on the full sky. It 
is common to treat this complication by a proper reduction of the potentially available 2.^-1-1 degrees of freedom at 
each £. 

We shall now see how these factors expand the confidence intervals of various estimators, and in particular the 
DJ^. Only by taking all these factors into account, wc will be in the position to test the 'null' (no gravitational 
waves) hypothesis (j33p . (j34p . and make predictions for the forthcoming experiments. 
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A. Effects of noises 



In general, the multipole coefficients a^„(o) observed in a CMB measurement consist of a signal convolved with the 
beam window fmiction, a^j(s)W£, and a collective noise contribution af„-^{n): 

For WMAP and Planck satellites, the window function Wi is close to 1 at ^ < 100 [H, SJ, so we set We = 1. 

As usual [2^ m, [13, we assume that the noise terms af^{n) and the signal terms af^{s) are statistically 
independent of each other. Also, any two members in the set {cif^{n)\X = T,E,B;i = 2,3,- • ■;m = —i,- ■ ■,£} 
are statistically independent, while each member af^{n) obeys a normal distribution with zero mean and standard 
deviation erf [n) , 

It is presumed that the noise power spectra Nf--^ are known from independent evaluations. Obviously, the TE noise 
power spectrum vanishes, Nf^ = 0, due to the absence of correlations between ajjy^{n) and af^{n). 



1. Probability density functions for multipole coefficients af^[o) and estimators Df'^ (o) 

Since and cif^{n) are statistically independent and obey zero-mean normal distributions, their sum afjy^{o) 

obeys a zero-mean normal distribution with the variance {erf (o))'^, which is the sum of the signal and noise variances, 

(af (o))2 = {af^{o)af^{o)) = {al,{s)af^{s)) + (afM^fM) = + A^f ^. 

The joint pdf f {ciJ„i{o) , af^^{o)) is a zero-mean bivariate normal distribution with the correlation coefficient Pi{o), 

^J{-UoH:,{o)){af^{o)af*{o)) ^ + NDiCf^ + N^) 

Comparing Eq. ((38|) with Eq. (fT9|) one finds that the noises make the output correlation coefficient pi{o) smaller than 
the no- noise correlation coefficient pi. 

The unbiased auto-correlation estimators in the presence of noise are defined by 



The expectation value and the standard deviation amount to 



(39) 



Comparing with Eq. ((23|) one finds that the noises make the standard deviation ADf-^{o) larger than the no-noise 
value ADf^. 

The pdf of the estimator Df^{o) with n degrees of freedom can be derived in a way similar to that in Sec. IIV Al 



n 



(T/(o))"/2-ie-^(°)/2 



2«/2r(n/2)(af (o))2 ' 



where V{o) = n{Df^{o) + A/^)/(af (o))2. 

Finally, we have to consider the unbiased cross-correlation estimator Dj^{o) which is of prime importance for us, 



1 ^ 

^^(0)= 3(3^^^) E (aL(o)af„:(o)+aF;;(o)af„(o)). 
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FIG. 6: The estimation of the WMAP5 noises Nf^ and TVf ^. 



The expectation value and the standard deviation amount to 

{Dj^{o))=pdo)aJ{o)af{o) 



'TE 



{Cji-+NjT){CfE+NEE) + (cJ'^y 
21+1 



Again, the noises make the standard deviation swell in comparison with the no-noise case ([28 
The pdf for the estimator Dj^{o) with n degree of freedom is similar to that in Sec. IIVB| 



Drio)\ 



X exp 



Piio) 



1 



1 



l-{p,io)yaJ{o)af{o) 



Kr 



7r(l-(p,(o))2)r(|) 

n\Drio)\ 
Xl-{p,{oW)aJ{o)af{o) 



(40) 



(41) 



To summarize, the pdf (j4ip . which includes noises, retains the functional form of the no-noise pdf (|3ip . but with the 
anticipated substitutions. Namely, Dj^ , aj, af , pg are being replaced with Dj^{o), aJ{o), af (o), pg{o), respectively. 



2. Numerical values of the noise power spectra 

The symbolic, so far, noises Nf^ have specific numerical values in concrete experiments. We use evaluations 
quoted in the WMAP and Planck literature (we adopt the experimental parameters but not the fundamental physics 
in these references) and apply them to the TE estimators in the WMAP5 best-fit model 

The evaluation of the 5-year WMAP noises, which include the instrumental noise, point sources noise and some 
systematic errors 'so], follows from Ref. [Ill, 0, |4^, IZtI. Starting with numerical values for ADj^ , ADj^ cited 
in [ill, and using Eq.(lO) from ^] and Eq.(l) from [47t. we derive the noise terms Nf^ , Nf-^ and plot them in 
Fig. [HI We use these graphs when we refer to the WMAP5 noises. The noise power spectra have a weak multipole 
dependence. For example, when i = 10, one has N^'^ = 1.46 x lO'^/xK^, Nf^ = 2.42 x lO'^^K^; while for £ = 100, 
one has Nj"^ = 1.36 x 10~^^K^, Nf^ = 2.59 x 10~^/zK^. (Actually, it seems to us that these numbers somewhat 
overestimate the true WMAP5 noises, as will be discussed in Sec. IVIII ) 
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FIG. 7: The correlation coefficient piip) (left panel) and the probability \og^Q{Dj^ {o) < 0) (right panel) in the WMAP5 best-fit 
model. In both panels the solid, dashed and dotted lines mark the no-noise, Planck and WMAP5 curves, respectively. 



For the Planck experiment, we use the expected instrumental noise in the 143GHz frequency channel. This channel 
is supposed to have low foreground levels and the smallest noise. The quoted noise power spectra have no ^-dependence 

TVj^ = 1.53 X 10"VK^ Nf^ = N^^ = 5.58 x lO^^K^- (42) 

It is demonstrated in Eq. (j38p that noises degrade the correlation coefficient pi{o). In the left panel of Fig. [7] we 
plot plio) as a function of ^ for WMAP5 and Planck experiments in comparison with the no- noise graph from Fig. [1] 
In the right panel of Fig. [7] we show P{Dj^{6) < 0) for WMAP5 and Planck experiments in comparison with the 
no- noise graph from Fig. [5] Clearly, large WMAP5 noises make the probability P{Dj^{o) < 0) dangerously close to 
0.5. 

In Fig.il we show the 68.3%, 95.4%, 99.7% confidence intervals of the Dj^{o) estimator for Planck and WMAP5 
experiments, starting from the no-noise case in the left panel (the same as the left panel in Fig. d]). Obviously, with 
increasing noise the confidence intervals become progressively wider. 



B. The effects of the cut sky 



Some additional stretching of confidence intervals should be attributed to the incomplete sky coverage. Although 
the shape of the sky cut is, in principle, important for accurate calculations (49l |. we adopt a simplified approach 
wherein the available number of degrees of freedom is reduced to n = (2£+ l)/sky, where /sky is the observed portion 
of the sky. We use the extra symbol c to denote pdf's in this cut sky approximation. Referring to Eqs. (|3T|). (|4T|) we 
write 

UDf^'io))^ fiDf^'io)) ,^ ^ . (43) 

n=( 2^-1-1) /sky 

The mean value of the XX' estimator calculated with fc{Df^ (o)) is the same as in the case of the full sky coverage, 
but the standard deviation increases, 

, , (44) 
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FIG. 8: The 68.3%, 95.4% and 99.7% confidence intervals in tlie WMAP5 best-fit model. Panels from left to right: no-noise, 
Planck, WMAP5. Note changes in scaling on the vertical axis. 



Comparing Eq. (|44|) with Eqs. ([39l [40|) we see that the standard deviation increases by the further factor l/-\//sky. 
In particular, the confidence intervals in Fig. [5] should be stretched in the vertical direction by this factor. 
The WMAP Collaboration [111 and the 'Planck Blue Book' [3 quote the cut sky factor 

/sky(WMAP) =0.85, /sky(Planck) - 0.65, (45) 

respectively. The former one leads to the 1 / -\/ /sky ~ 1-08 increase in the uncertainty of the estimators, and the latter 
one leads to 1/ \/ fsky ~ 1.24. 

In our further analysis we use these stretched confidence intervals, but it may turn out that the quoted noises and 
/sky are overly pessimistic, especially for Planck, so the real level of uncertainty may prove to be smaller. 



VII. TESTS OF THE NULL HYPOTHESIS: WMAP5 BEST-FIT MODEL WITH NO GRAVITATIONAL 
WAVES 

The developed theory can now be confronted with real and simulated data. In Fig. [51 based on information from 
Ref. [l,|4l|, we show the WMAP5 unbinned TE data points on the background of the derived in Sec. lVIi r_E confidence 
intervals. Despite the messiness of the picture, the visual impression is such that the curve of the 'center of gravity' 
of the data points lies somewhat below the black solid curve expected of density perturbations alone. This tendency 
is of course required by the presence of gravitational waves, as we discussed in Sec. [V] The visual impression should 
be quantified by means of the 'null hypothesis testing'. The null hypothesis under investigation, denoted iJo, is the 
WMAP5 best-fit model (|33p . ([M)) with no gravitational waves. We probe Hq with the mean value (M) and the 
variance (K) tests. 



A. The mean value test 

We introduce the statistic M 

Dr{o)-cr 



To test observations, one takes Dj^{o) as the actually observed values of the estimator. The expectation amount 
Cf^ and the standard deviation ADj^{o), which weighs the scatter of the data points in Eq. (H^ . are calculated 
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FIG. 9: The WMAP5 unbinned TE data points on the background of confidence intervals which include WMAP noises from 
Fig. [6] and /sky = 0.85. The black solid line is the expectation value according to the Ho hypothesis. 



from Eq. (|44ll assuming the validity of the hypothesis Hq. The negative (positive) value of M quantifies the propensity 
of the data points to lie below (above) the theoretical prediction. For unbinned data, the summation in Eq. (|46p is 
over £ = 2,3,- • -,100. 

In confirmation of the visual impression, the value of M on the actually observed TE data points turns out to be 
negative and equal to 

TO = -8.69. (47) 

To assess the case for rejection of Hq one needs to know the pdf for the variable M . Using pdf ([^5]) we generate 10''' 
random samples of the estimator Dj^{o), where each sample includes all £ — 2,3, ■ ■ ■, 100. We then calculate M for 
each sample. The distribution of these M (effectively, the pdf for the variable M) is shown in Fig. [TOl 
The significance level a of the test, or in other words the probability of incorrect rejection of Hq, is 

P( \M\ > ItoI 



Ho 

Calculating a from the pdf in Fig. 1101 and with m from Eq. (j47p . we find 

a = 40.63%. (48) 

Although the case for rejection of the hypothesis Hq is weak, within la, Eqs. (j47l [48|) serve as an indication that one 
may try to find a better hypothesis. 

To check the stability of this conclusion, we expanded the M test. First, we carried out the same procedure with 
the 3-year WMAP data and found similar results. Second, we assumed the full sky coverage /sky = 1.0 with WMAP5 
unbinned data and arrived at a slightly stronger evidence for rejection of Hq, to = —9.43, a = 32.91%. Third, we 
applied the M test to the WMAP5 data binned at thirteen multipoles i = 4J_0JJ, 22, 27, 33, 40, 48, 56, 65, 76, 87, 98 
[jj EH , as shown in Fig. [TT] We made simplifying theoretical assumptions [H, 0, [sl] according to which the estimator 
at every binned multipole obeys a Gaussian distribution with the mean value Cj^ taken at £ representing the bin, 
and with the standard deviation properly averaged over the bin and reduced by the size of the bin. Binned multipoles 
are supposed to be statistically independent. In Fig.[TTl the binned WMAP5 TE data {£+ l)Dj^(o)/27r are shown by 
(blue) dots, and the la intervals, which include WMAP5 noises from Fig. [6] and /sky = 0.85, are shown by the largest 
(blue) bars. Following a somewhat incorrect practice, we transferred the la confidence intervals from the theoretical 
pdf to the data points themselves. The size of the bins is indicated by the width of the (yellow) rectangular regions. 
The black solid line is the expectation value according to Hq hypothesis. Under the simplifying assumptions made. 
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FIG. 10: The outcome values ol the test quantity M in simulated samples. In the actually observed realization, m = —8.69. 
TABLE L Significance level a in the mean value tests based on WMAP5 data. 



data 


noises 


a (cut sky) 


a (full sky) 


unbinned data 
binned data 
binned data 
binned data 


WMAP5 noises 
WMAP5 noises 
Planck noises 
no noise 


40.63% 
30.21% 
1.39 X IQ-^ 


32.91% 
26.52% 
7.08 X 10"** 
5.19 X 10"" 



the quantity M satisfies a zero-mean Gaussian distribution with variance n — 13, determined by the number n = 13 
of the observed (binned) data points. The analogue of Eq. (|46|) . with summation over thirteen binned muhipoles, 
resulted in to = —3.72, a = 30.21%. The additional assumption /sky = 1 reduces a to a = 26.52%. 

The crucial importance of efforts to reduce noises in future experiments is illustrated by the foUowing exercise. We 
assumed that the Planck mission would observe exactly the same TE data points as WMAP5 did, but which should 
now be analyzed against the Planck's smaUer noises. We first calculated the confidence intervals of the unbinned 
data using the Planck noises (|42|) and/sky = 0.65, and then we have binned the data using the same procedure that 
was used in the WMAP analysis [1, |4lj. The resulting la intervals are shown by the (red) intermediate-size bars in 
Fig. [m With these smaller uncertainties, the M test gives to = —15.67, a = 1.39 x 10"^. The further assumption 
/sky = 1 makes a even smaller, a — 7.08 x 10"^. In other words, in these circumstances, the null hypothesis would 
be rejected at a higher than 4tT level (a ~ 7 x 10"^). 

Finally, we have looked at what would happen in the ideal, no-noise, case, that is, when the standard deviation 
is at its minimum level AI^J^ and /sky = 1. The la intervals are shown by the height of the (yellow) rectangular 
regions in Fig. [TT] The M test results in m = -27.17, a = 5.19 x 10"^''. Thus, in the ideal case, the WMAP5 TE 
observations would have rejected the WMAP5 Hq hypothesis at a higher than 7a level. The results of all M tests are 
summarized in Table HI 



B. The variance test 



The scatter of the observed data points is probed by the K statistic 
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FIG. 11: In this figure, based on 5-year WMAP observations 0, |4l[, we show binned TE data points {£ + 1)DJ'^ {o)/2n 
surrounded by various la uncertainties: WMAP5 (blue, largest), Planck (red, intermediate), no-noise (yellow, smallest). 



The significance level a of the K test is defined by 



P K > k 



He 

First, we apply the K test to the WMAP5 TE unbinned data, shown on Fig. (5] The pdf of the test statistic K, buih 
by the simulation method described in the previous subsection, is shown in Fig. [T^l The value of the K statistic and 
the corresponding significance level are as follows k = 85.38, a — 82.63%. Under the assumption /g^y = 1, a reduces 
to a = 44.38%. Thus, with WMAP5 noises, the variance test does not provide enough evidence for rejection of the 
null hypothesis. Nevertheless, the somewhat lower than expected value of k points toward a possible overestimation 
of the WMAP5 noises. 

We have also applied the K test to the binned 5-year WMAP data. We have made the same simplifying Gaussian 
assumptions that were discussed in Sec. IVII Al Under these assumptions, the quantity K obeys a Xn distribution 
with n = 13 degrees of freedom [5^. The peak of this function is at X w n. Using this pdf, we find k = 11.48 and 
a = 57.06%. The full sky assumption reduces a to a = 40.96%. Again, with WMAP5 noises, the case for rejection 
of Hq is weak. 

The value of a dramatically decreases when the WMAP5 data are accompanied by the Planck noises. In this case, 
k « 211 and a ~ 10~^^. The Hq hypothesis would have been rejected with a huge margin. The a is even smaller in 
the no-noise case with /sky — 1- The summary of K tests is shown in Table HIl 

It is important to note that, with the Gaussian approximations made for the pdf of the estimator Dj^{o), the 
variance test is equivalent to the test, where the calculations are being done directly with the weighted values of the 
pdf at individual data points: 



is: = -2 In 



TE\ 



21nn 



fiprio)) 



'TE\ 



(49) 



where f{Dj^{o)) is the value of the pdf at the observed data point Dj^{o), and f{Cj^) is the value of the pdf at 
the expected point Dj^{o) = Cj^ . 



VIII. LIKELIHOOD ANALYSIS OF THE QUADRUPOLE RATIO R 



The existence of relic gravitational waves is a necessity dictated by general relativity and quantum mechanics, 
and not an inflationary 'bonus' which, as follows from inflationary theory, most likely should not be awarded, see 
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FIG. 12: The distribution of the test quantity K from the simulated samples. In the observed realization, k = 85.38. 
TABLE II: The significance level a oi K tests of WMAP5 data 



data 


noises 


a (cut sky) 


a (full sky) 


unbinned data 
binned data 
binned data 
binned data 


WMAP noises 
WMAP noises 
Planck noises 
no noise 


82.63% 
57.06% 
« 10-3^ 


44.38% 
40.96% 

< 10"^'' 

< 10"^'^ 



Sec. im We do not use our theoretical position as a technical 'prior' in data analysis, but we certainly are willing to 
reconsider the WMAP5 TE data from this point of view. We also intend to make predictions for future experiments. 
The immediate goal is to evaluate the quadrupole parameter i?, Eq. p6p . from the existing observations. The mild 
indications in the WMAP5 TE data favoring the rejection of the Hq hypothesis (i? — 0) have been discussed above. 



A. Expanding the parameter space 



The inclusion of gravitational waves with two new free parameters and rtj (see Eq. ^) would require, strictly 
speaking, a new full-scale likelihood analysis of all the data. Certainly, the addition of gravitational waves must be 
consistent with all the measured correlation functions and upper limits, Eq. (|15p . In what follows, we simplify this 
approach, while retaining the main points of what we want to demonstrate. First, we keep intact the background 
cosmological parameters, as quoted in Eq. ([33|) . Second, we bound the four perturbation parameters B^, ris, -B^, nt 
by three constraints, thus leaving independent only one parameter. We choose this independent parameter to be the 
quadrupole ratio R. The three mentioned constraints are as follows: 

nt = ns - 1, + l)CUio/27r = MQ^iK^, ns{R) = n^O) + 0.35i? - O.OTi?^ (50) 

The first constraint is purely theoretical, it reflects the origin of cosmological perturbations [i^ . The second constraint 
plays the role of an overall normalization, it requires the perturbation parameters to satisfy the demand that the 
joint, dp plus gw, temperature anisotropy at ^ = 10 were fixed at the well-measured level. The third constraint, 
where ns(0) = 0.960 [J|, is empirical. The function ns{R) was designed to be such that the best studied correlation 
function Cf^ were as close as possible to the actually observed data for a wide range of possible i? > 0. In contrast to 
other analyses, we are not using the inflationary 'consistency relation' nt — ~r/8, which affects the results regarding 
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FIG. 13: The TT (left) and EE (right) power spectra in models with varying amounts of gravitational waves. The dashed 
(red) line is for the Hq hypothesis, R = 0. The solid (blue) line shows the model with R — 0.24, and the dash-dotted (green) 
line is for the model with R = 0.3. 

gravitational waves (see comments in j5| on the self-contradictory nature of the data analysis based on inflationary 
formulas) . 

In Fig. [T3] we iUustrate, for selected values of R, the consistency of our class of models with the TT and EE 
observations. As usual, some corrections at very high £'s could be achieved, if desired, by the assumption of the 
'running' spectral index Ug, but we did not resort to this possibihty. It is also important to emphasize that the 
BB spectrum in our class of models is well below the existing upper limits [43]. Indeed, the WMAP limit is £{i + 
l)C£^_g/27r < 0.15^J.K'^ (95%C.L.). This corresponds to r < 20 0], which is roughly R < 10, wher eas we will be 
mostly dealing with R < 1. 

B. Likelihood function for R from the WMAP TE data 

The likelihood function for a parameter is a pdf for a random variable, where the parameter is considered unknown, 
whereas the value of the random variable is considered known from experiment (see for example 29. l49l. ISTj). Since 
the variables Dj^{o) with differing £'s are independent, their joint pdf is a product of individual pdf's fc{Dj^{o)) 
over £. Therefore, the likelihood function for R is 

C{R) = C\{UDr{o)), (51) 

I 

where fc{Dj^{o)) is given by Eqs. ([43l [4T|) . the quantities Dj^{o) are taken from observations, and C is a normal- 
ization constant to be fixed later. 

Moving by small steps Ai? = 0.01 from i? = we build C{R) numerically. The Dj^{o) participating in Eq. ([?T|) 
are the unbinned data shown in Fig. [HI the other quantities are determined by Eqs. (|331 ISO)) . Fig. [S] and /sky — 0.85. 
The result of this calculation is shown by a solid line in Fig. [Ml 

The maximum of C{R) is at i? = 0.240. The constant C was chosen to make C{R) = 1 at maximum, as illustrated 
in Fig.[T31 Taking into account R — 0.240 and Eq. ([SO)) , the full set of our best-fit perturbation parameters is defined 

by 

e{£ + l)CUio/27r = 840AiK2, R = 0.240, = 1.040, nt = 0.040. (52) 

This model has a blue(ish) primordial power spectrum and accommodates a significant amount of relic gravitational 
waves. The TE spectrum for this model is shown by a solid line in Fig. [151 The main difference with the TE 
spectrum of the WMAP5 best-fit model (dashed line) is at multipoles £ < 50. This could be anticipated, as the role 
of gravitational waves is essential only at ^ < 100. The TT and EE spectra for our model (i? = 0.240) are plotted in 
Fig. [HI 

The difference with WMAP conclusions is indicative but not significant, owing to large WMAP noises. Indeed, the 
confidence interval around the maximum likelihood value R = 0.240 is broad. Using the techniques of Appendix [^ 
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FIG. 14: The likelihood function for R. The solid line uses the unbinned 5-year WMAP data, whereas the dashed line is for 
the unbinned 3-year WMAP data. 



we find that the 68.3% confidence interval is represented by i? = 0.240^0 225- (Lower boundary cannot be less than 
—0.240 as the quantity R is non-negative.) Thus, the value i? = is only barely outside the 68.3% confidence interval. 
This result is consistent with the null hypothesis testing in Sec. IVIII 

We have also calculated the likelihood function C{R) for 3-year WMAP TE data. The background cosmological 
paxameters are taken as the WMAP3 ACDM model: fl^/i^ = 0.1265, fJfc/i^ = 0.0223, h = 0.735, r^ejon = 0.088 
jsil. In line with the WMAP3 findings we set ns{0) = 0.951 [El] and use Eq. ((SO]) to relate Us and R. The resulting 
likelihood function is shown in Fig. [T3] by a dashed line. The maximum likelihood value, along with the associated 
68% confidence interval, is i? = 0.149t°;^^^. The best-fit spectral indices are Ug = 1.002 and nt — 0.002. Again, the 
i? = is within the 68.3% confidence interval. Nevertheless, it is worth noting that in transition from WMAP3 to 
WMAP5 data, the peak value of R increased from R = 0.149 to R = 0.240. 



C. Likelihood analysis of simulated TE data 



Being guided by the experience that in astrophysics very often 'today's indication is tomorrow's discovery', we treat 
our best-fit model (|52p as a benchmark model to be confirmed by Planck and other sensitive missions. Using the 
parameters of this model we simulate the TE data for Planck, taking into account the expected noises and /sky We 
also vary the perturbation parameters, noises and sky coverage in order to assess what will be happening with nearby 
models. 

The procedure for Planck is as follows. First, we choose the input value i?i and determine other perturbation 
parameters from Eq. ([50)1 . The background parameters are fixed by Eq. ((33|) . Second, we adopt Planck parameters, 
Eqs. (1421 H5|) . and randomly generate 10 sets of data points {Dj^{o)\i = 2,3, • • •, 100}, where each set includes all 
indicated ^'s. We use for this purpose the underlying pdf/c(£'f^(o)), Eqs. (I13[1I|). Obviously, even if all parameters 
are fixed, the outcome data points Dj^{o) are supposed to be random due to the randomness of signal and noises. 
Third, for every set of simulated data we build, in a manner described in Sec. IVIII Bl the likelihood functions C{R) 
and explore their properties. 

Let us start from the benchmark value i?i = 0.24. The 10 likelihood functions built from 10 random realizations are 
shown in Fig. [1^] (left plot). As expected, the maxima of the likelihood functions concentrate around the input value 
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FIG. 15: The solid line is our best fit to the WMAP5 binned TE data, in comparison with the Ho hypothesis (dashed line). 
The depicted error bars are taken from the R = model. In the insert, we plot the same graph for unbinned data. 



Ri ~ 0.240, and the widths of the hkclihood functions are approximately equal to each other. The arithmetical mean 
of the maxima is Rp = 0.242, and the arithmetical mean width of the 68.3% confidence intervals, i.e. la interval, is 
AR = 0.074. The realization whose maximum is closest to the input value Ri = 0.240 (in a sense, a 'typical' likelihood 
function) is shown by a dark line. 

On the ground of these simulations we conclude that the Planck satellite will be able to see a relatively strong 
evidence for relic gravitational waves, if the Ri — 0.240 model is correct. In other words, Planck will be able to 
exclude the value i? = at more than 3cr confidence level. We have also done simulations for a more optimistic cut 
sky factor /sky = 0.85. In this case, AR — 0.062, thus allowing to exclude the R — assumption at 4a level. 

The best achievable conditions for detecting the signature of relic gravitational waves are reached in the idealized 
case of no noise and full sky coverage. With the input value Ri = 0.240, the 10 likelihood functions for TE experiment 
in ideal conditions are shown in Fig. [16] (right plot). The mean value of the maxima is Rp = 0.228, and the mean 
width of the 68.3% confidence intervals is AR — 0.032. Thus, the ideal TE experiment would be able to detect 
gravitational waves and exclude i? = at more than 7a confidence level. 

We have also evaluated the R detectable at a 2a level by Planck or ideal experiment. In the left panel of Fig. [T7] 
we show 10 likelihood functions for the input Ri = 0.110 and Planck parameters. The mean value of the maxima is 
at Rp = 0.108, and the mean width of the la confidence intervals is AR = 0.067. Thus, using the TE method, the 
Planck satellite would be able to detect R = 0.110 at nearly 2a level. The left panel of Fig. [T71 shows the results of 
solving the same problem in the ideal case. The input R for the 10 likelihoods is Ri = 0.050. The mean value of the 
peaks is Rp = 0.051, and the mean width of the la confidence intervals is AR = 0.025. Thus, an ideal experiment 
would detect R — 0.05 at exactly 2a level. Even more possibilities are considered below in the context of comparison 
of TE and BB methods. 



IX. COMPARISON OF TE AND BB METHODS IN SEARCH FOR RELIC GRAVITATIONAL WAVES 
AND PROSPECTS FOR THE PLANCK MISSION 

The advantage of the usually discussed B-mode searches for primordial gravitational waves is that in the BB 
correlation function there is no competing contribution from density perturbations. However, the BB signal is 
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FIG. 16: The likelihood functions for R with the input value _Ri = 0.240. The Planck conditions are used in the left panel, the 
ideal conditions - in the right panel. The 'typical' likelihood functions are shown by the dark lines. 




FIG. 18: Likelihood functions for R from simulated BB data with Planck parameters Nf^ — 5.58 x 10"^fiK^, /sky = 0.65. 
The input values in the panels from left to right are Ri = 0.24. 0.07, 0.04. The dark solid lines show the 'typical' likelihood 
functions. 

inherently weak, it is about 50 times smaller than the g.w. contribution to the TE correlation. It is true that the 
nominal BB noise is also smaller than the TE noise, so that the nominal 'signal to noise' ratio may be slightly in 
favor of the BB method. But one must take into account the large systematic effects that the i?-mode is prone to. 
They include the various contaminations [H, [s^, [g^, [6l[ , the disadvantages of dealing with an auto-correlation {BB) 
rather than with a cross-correlation {TE) function, the actual weakness of the signal which can lead to the danger of 
unexpected leakages and couplings, etc. It is only the experimenters who can properly assess these complications. 

When only the evaluation of the Planck's instrumental BB noise, Eq. (|42|). is used in calculations, we call it 
an 'optimistic' BB case. As for the additional complications mentioned above, we try to treat them as 'effective' 
noises in the framework of our formalism. We give our evaluation of these 'effective' noises and increase the nominal 
instrumental BB noise to the level which we call the 'realistic' BB case. We think that it is only this 'realistic' BB 
case that can be taken as a fair comparison with the TE method. The comparison of TE and BB methods was also 
considered in a recent paper js^] (see also earlier papers [62l. [63j). 

A. The likelihood function for R in the BB method 

We start from the 'optimistic' BB case where the BB noise is taken as Nf^ = 5.58 x 10~'*/.tK^, see Eq. p2)) . We 
generate the mock BB data using the pdf fc{Df^{o)) from Eq. Eq. ^ (setting XX = BB), and Eq. ((i5|) . 

The procedure is exactly the same as was described in Sec. IVIII Cl in the context of the TE mock data. The 10 sets 
of simulated data for every input Ri are then used in the likelihood function 

C{R)^CY[UDf^{o)), 

i 

where C is a normalization constant. 

First, we discuss the benchmark model Ri = 0.240. The 10 likelihood functions are shown in Fig. [18] (left panel). 
The mean value of the maxima is Rp — 0.237. The mean width of the 68.3% confidence intervals is AR = 0.050. This 
is a bit smaller than the mean width in the TE likelihood functions. For illustration, we considered other examples 
as well. The 10 likelihood functions with the input value Ri = 0.070 are plotted in the middle panel of Fig. [THl The 
mean values are Rp = 0.070 and AR = 0.024. This shows that R = 0.070 could be detected by the BB method 
at nearly 3a level. Finally, the likelihood functions with the input value Ri = 0.040 are shown in the right plot of 
Fig. [THl The mean values are Rp = 0.038, AR = 0.018. This R is at somewhat better than 2a level of detection. 
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TABLE III: Comparison of the TE and 'optimistic' BB methods in terms of the signal to noise ratio S/N . 





0.30 


0.24 


0.20 


0.15 


0.11 


0.07 


0.04 


TE: Rp 


0.28 


0.24 


0.19 


0.13 


0.11 


0.07 




TE: AR 


0.08 


0.07 


0.07 


0.07 


0.07 


0.07 




TE: S/N 


3.90 


3.24 


2.82 


2.21 


1.64 


1.06 




BB: Rp 


0.31 


0.24 


0.20 


0.16 


0.11 


0.07 


0.04 


BB: AR 


0.06 


0.05 


0.04 


0.04 


0.03 


0.02 


0.02 


BB: S/N 


5.36 


4.80 


4.55 


3.95 


3.44 


2.92 


2.22 



TABLE IV: Signal to noise ratio S/N in the 'realistic' BB case. 



Ri 


0.30 


0.24 


0.20 


0.15 


0.11 


0.07 


BB: Rp 


0.29 


0.26 


0.20 


0.15 


0.10 


0.06 


BB: 'AR 


0.13 


0.11 


0.09 


0.07 


0.06 


0.05 


BB: S/N 


2.40 


2.24 


2.30 


2.06 


1.72 


1.32 



Introducing the quantity 

S i?i 

as the 'signal to noise' measure of detecting gravitational waves, we show in Table Hill the comparative performances 
of TE and 'optimistic' BB methods in terms of S/N for a range of input values Ri. The Table also shows the mean 
values -Rp of the maxima of the likelihood functions and the mean widths AR of the 68.3% confidence intervals. 
(The number of displayed decimal digits is an artefact of numerical calculations, and not a demonstration of our 
responsibility for the accuracy of results.) 

Moving to the discussion of the 'realistic' BB case, we think that all the non-instrumental 'effective' noises can 
increase the level of the noise variables cif^{n), Sec. IVIl by at least a factor of 2. This means that the noise term 
Nf^ participating in our analysis, should be raised from its level in Eq. by at least a factor of 4. In other words, 
we are working with the 'realistic' noise term N^^ = 2.24 x lO^^/iK'^. Repeating all the calculations with this noise 
term, we come up with the 'realistic' evaluation of the Planck's BB performance, as shown in Table [TVl 

The numerical results for TE and BB methods are also presented as individual S/N points in Fig. [121 The solid 
lines in this figure show the analytical approximations which we are set to discuss. 



B. Understanding the signal to noise ratio S/N in TE and BB methods 

Our numerical calculations are supported by simple analytical expressions for S/N in TE and BB methods. Let 
XX' denote either TE or BB. With Gaussian approximation for individual pdfs, Eq. ((37l) . and given the statistical 
independence of differing £'s, the likelihood function C{R) reads 



£{R) = n 



1 



2ttAD 



XX' 



exp 



1 / D 

2 



XX' 



-C. 



XX' 



AD 



XX' 



The i?-dependence of C{R) is coming from the i?-dependence of Cf^ and ADf-^ . 

The estimate of the width AR of C{R) with the input (maximum likelihood) value i? = _Ri is given by 



AR = 



/ d^\nC{R) 




\ dR^ 





(53) 



where the angle brackets denote the average over the joint pdf for Df-^ . The second derivative participating in 
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we arrive at an intermediate expression 
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(54) 



Since the noises are i?-independent and dADf-^ /dR is smaller than dCf^ /dR by a factor oc l/-\/ {2£ + l)/sky, one 
can neglect the second term in Eq. (|54p bringing it to 



AR = 
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e=2 



da 
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AD 



XX' 



dR 



(55) 



The general expression (1551) is a well-known result [13] applied to CMB subjects on several occasions (see for example 
[i^ll^). Eq. ([55]) allows one to write the following analytical estimate for S/N: 



S/N=^ = 
' AR 



\ 



1=2 



1 acf^' 

AD^^' d\nR 



(56) 



The further simplifications of Eq. (|56[) rely on further assumptions. It is reasonable to assume that for sufficiently 
small R one can use the Taylor expansion 
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Then, for XX' = TE one derives 



S/N = 



\ 1=2 \ 



-R=0 



AD 



TE 



and for XX' = BB one derives 



S/N 



\ 



y- i CriB) 

i=2 
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BB 



(58) 



(59) 



These expressions could be expected on physical grounds. In the TE case, the quantity Cj^l „ „ is not zero, and 



\R=o 



whereas the i?-noise is given by the 
is zero (no gravitational waves - no BB 



it is entirely determined by density perturbations. Therefore, it is reasonable that the i?-signal, allowing to detect 
gravitational waves, is given (at each £) by the difference Cj^{R) — Cf^ 
total standard deviation ADj^ . In the BB case, the quantity Cf^ 

correlation), and therefore it is reasonable that the i?-signal is given by Cf^{R), whereas the i?- noise is given by 
ADfB_ 

As one could anticipate, the linear approximation ()57|) is not quite good if R is not small. Some deviations from 
exact numerical results could be expected. However, in Fig. 1191 we plot the TE curve without any corrections, that 
is, as it follows from the analytical formula ([55)1 . As for the BB case, the noticable deviations exist, so a better fit 
to the exact numerical results requires a small correction to Eq. (|59|) . The actually plotted BB curves in Fig. [19] are 
given by the corrected analytical formula 



S/N 



\ 



1=2 



Cf^iR) 



(l + 0.5i?) AD 



BB 



(60) 



The analytical results shown in Fig. [11] are based on imax = 100. 

It is natural that the total [S/NY consists of individual contributions {S/N)1 at each i. Both, Eq. ((SS]) and 
Eq. ([60]) . have the structure 



{S/Nf^Y.{S/N)l 



1=2 
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FIG. 20: The BB and TE gravitational wave signals (power spectra) in comparison with noises (standard deviations) in the 
benchmark model R = 0.24. The signals are in black, the noises are in red. 



where the individual signal to noise ratios {S/N)£ are given by 

fC/l\T\TE _ '^I'^jR) - Cf^jR = 0) fQ/isT\BB _ (1 + 0-5i?)^^Cf^(i?) , . 



One can say that the numerators in Eq. (j61[) are signals depending on R, while the denominators are noises. Signals 
and noises are shown, as functions of £, in Fig.[20]for our benchmark model with R = 0.24. The lower solid black line 
for the BB power spectrum includes the correcting factor (1 + 0.5i?)~^, whereas the upper solid black line does not. 

Fig. [20] provides a graphical illustration of what the signal to noise ratio S/N in Fig. [10] consists of and how it is 
being accumulated from the individual, predominantly lower than the noise, contributions at each i. It is seen from 
Fig. [20] that the success of the BB method strongly relies on the very low multipoles, £ < 10, and hence on the era 
of reionization with the optical depth t from Eq. (jSH)) . Most of S/N in the BB method comes from these multipoles 
Smm. Specifically, in the R = 0.24 model we find 

The cosmic reionization remains to be fully understood and quantified (67j . If it happens that the optical depth is 
actually smaller than the currently accepted Trdon ~ 0.08, the sensitivity of BB method to relic gravitational waves 
will be reduced. On the other hand, if the BB detection does take place, it will be an argument for the currently 
accepted or higher r. 

In contrast to BB, the TE method does not rely on very low multipoles. In fact, the most of S/N comes from the 
multipoles 30 < £ < 70. Specifically, in the R = 0.24 model we have 



Y.t2{S/N)l 

El%{s/N)l 



0.15. 



One can say that that the TE correlation function directly probes the relic gravitational waves at recombination. 

It emerges from this comparative analysis and Fig. [TO] that both, TE and BB, methods need to be used, with the 
former one being more likely to bring the first identification of relic gravitational waves in the Planck data. 
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X. CONCLUSIONS 



A correct theoretical picture of the physical phenomenon is important for a proper setting of search experiment, 
as well as for the analysis of data and its interpretation. The relic gravitational waves 'must be out there', and 
we take this understanding as our starting point. Having further developed statistical theory of CMB anisotropics, 
we concentrated on the TE correlation function. We have shown that the WMAP5 TE data contain a hint on the 
presence of gravitational waves in the data. Our best-fit model includes a substantial amount of relic gravitational 
waves, R = 0.24. In simple terms, this means that 20% of the temperature quadrupole is caused by gravitational 
waves, and 80% by density perturbations. Because of large WMAP5 noises, the confidence of this conclusion is not 
high. The maximum likelihood result includes the WMAP's best-fit model with no gravitational waves, R — 0, almost 
within la interval. 

We projected our conclusion R = 0.24 on the forthcoming Planck mission. We numerically simulated Planck data 
and compared the TE and BB channels of observations. In the BB channel, much will depend on contamination by 
systematic effects. We distinguish the 'optimistic' BB case, when only the advertised instrumental noises are taken 
into account, and the 'realistic' BB case, when the effective noises are increased as a reflection of many possible 
residual effects. It is shown that the TE method will see the relic gravitational waves at a better than 3ct level, 
whereas the 'realistic' BB method will see them at a better than 2a level. 

Although most of our discussion focused on the Planck experiment, some balloon-borne [H, [6^ and ground-based 
experiments [lO, [lH [ll, currently in preparation, wih be sensitive to part of the lower-£ spectrum of CMB. 
The contributions of various intervals of i to the total signal to noise ratio S/N can be read off from the evaluations 
developed here. These experiments will be complementary to the Planck mission and can provide a healthy competition 
in the race for discovery of relic graviational waves. 
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APPENDIX A: PROBABILITY DENSITY FUNCTION FOR v 



In Sec nil B"3l it was shown that the joint pdf's /(-^ 



V2 ' V2 



(where to > 1) are bivariate zero- mean 



normal distributions with standard deviations ^ and correlation coefficient p^. The indices £, m, c are considered 
fixed, so the two variables are distinguished by the labels T, E. From such a bivariate distribution one can derive a 



pdf for the product variable 



T(c) E(c) 



(see Eq. 



For each w^^, the corresponding pdf is given by [74 



^af^ 



exp 



(1 - P^"" 



,)a^ cr. 



(Al) 



where Kq is the zero-order modified Bessel function. 

Let us analyze the pdf (|Aip in more detail. Firstly, Eq. (|Aip shows that the pdf's with positive and negative 
are related by 



(A2) 



If = (no correlations), the pdf is symmetric with respect to the axis v'f^^ = 0. For any p^, the pdf is divergent at 

^'fm ~ because the modified Bessel function KQ{a) is divergent at a = 0. 

In Fig. [21] we show the pdf's for different values of positive p^. When pi increases, the distribution shifts to the 
right refiecting the fact that the mean value (f^^) is proportional to pi. In the limiting case pi ^ 1 the pdf develops 



a step, making f{vf^) 

of negative, or positive, values of are given by 



for < 0. In general, using (6.621.3) from Ref. [35], one can show that the probabilities 



<0^ 



^(l-p2)l/2 



F 



1 3 p£ 

'2'2'p, 



(A3) 



34 



FIG. 21: Probability density lunctions /(«^^) for various positive pis. The solid, dashed and dotted lines depict pi = 
1.0, 0.4, and 0.0. 



P 



i:^ > o) 



^ (1 - p2)V2 



F 



T 1 ^ PL 



(A4) 



where F is the hyper geometric-iunction. These probabilities depend only on pi, but not on aj and af. 

Similarly to what has been done in [t^, one can show, using (9.131.1), (9.121.7), (1.624.9), (1.623.2) from [SE 



that -P(w£^ < 0) + -P(w^^ > 0) = 1, which confirms the correct normalization of /(f)^). As expected, when 



Piv^/J, < 0) = and P{v';^ > 0) 



1. And when pe = 0, P(v^;2 < 0) = P(w^;^ > 0) = 0.5. 



The probability density function (jAip allows one to define the confidence intervals of finding specific values of v^'^^ 



hn ■ 



As an example, we start from the 68.3% confidence interval. We denote by (w^^)[/ and ('I'^^j)l the upper and lower 



boundaries of the shortest interval, i.e. the interval having the minimum difference (w^^);/ 



(^£m)i' '^uch that (75 



fUm) d v'.^l - 0.683. 



The surface area under the pdf curve between (w^^)^/ and (w£^)l gives the probability (68.3%) of finding in this 
interval. This is a sensible measure [t^ of uncertainty of our estimator. For the pdf's that we are dealing with, the 
68.3% confidence interval always includes points of the maximum and the mean value of the pdf. 

In Fig. [211 we plot the 68.3%, 95.4%, 99.7% confidence intervals for different 's. When pi — 0, the confidence 
intervals are symmetric with respect to the mean value (zero) of w^^^j. When pf increases, the intervals move toward 



larger values w^^. In the limit 



1, the lower boundary (u^,^)l approaches for all levels of confidence. For any 



confidence level, the confidence interval is not symmetric with respect to the mean value of vl' . It can also be seen 



that, unless pi 
non-negative. 



(c) 

1, there is a significant probability of finding v^J^ < 0, despite the fact that the mean value is strictly 
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Correlation coefficient: 



FIG. 22: The 68.3%, 95.4% and 99.7% confidence intervals of tlie variable «^„/crJcrf for different pe. The solid (black) line 
denotes the mean value of v'^'^/ajaf . 



APPENDIX B: SOME PROPERTIES OF f{D] 



It is insightful to discuss properties of f{Df^) in comparison with /(w^^) in Eq. (|A1|) . Surely, both functions go 
to zero when the argument — > zLoo. But the behaviors at zero argument are different, f(Dj^) does not diverge at 
DJ^ = if n > 1 . Although the modified Bessel functions do diverge at zero argument jS^] , 



{v - 1)! (a\-'^ 



— ) , where {v>l, a ^ 0), 



the first factor in the right hand side of Eq. (PT]) compensates for this divergency, so that the pdf (I3ip remains finite. 

In the hmit ^ 1, keeping other parameters fixed and using the asymptotic formula 35] Ki,{a) ^ \/^' 
a ^ oo, one obtains 



With V = ^ p . 



2"/2r(n/2)aJcr^ 



As expected, the asymptotic distribution is a distribution with n degrees of freedom. 
Similarly to Eq. (jA2|) one has 



The probability of negative, or positive, values of Dj^ is given by 



p{Dr<o) 
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FIG. 23: PiDf^ < 0) as a function of e for pt = 0.0, 0.2, 0.4, 0.8, 1.0. In the upper left panel the vertical axis is PiDf^ < 0), 
whereas in other three panels it is logj^Q P{DJ^ < 0). 



or 



Pinr > 0) 



4-00 



r{n) 



r(f + i)r(f 
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n P£ + 1 



1,1-2,1 



-1 
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41) when 



These probabilities depend on pi but not on cr^ .crf . As it should be, Eqs. (|B11 IB2p reduce to Eqs. 
n = 1. 

We have checked numerically, for a set of ^'s and p^'s, the normalization condition P{Dj^ < 0) + P{Dj^ > 0) = 1. 
A useful equality relating the probabilities for positive and negative piS is 



= 1 - p{dJ^ < 0) 



A similar equality is valid for probabilities P{Dj > 0) | . When = 1 , we find, using formula — ^,1 + ^^,0] — 

1, that P{DJ^ < 0) = for any n. On the other hand, when pi = 0, we find, using formula F [l, 1 — ^, 1 + ^, — l] = 



that PiDj ^ < 0) = P{Dl ^ > 0) = i for any n 



2 r(^) ' ^ "-^e - - - \-^£ ^ - 2 
In Fig. [23l we plot the probability P{Dj^ < 0) in the interval 2 < £ < 70 for selected values of p^, pi — 
0.0,0.2,0.4,0.8 and 1.0. It is seen that for a fixed pi the value of P{Dj^ < 0) rapidly decreases as £ becomes larger, 
excepting the limiting cases p£ = or 1. For a fixed £ and growing pi, the P{Dj^ < 0) changes from 0.5 to 1. 
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